Una serie de tiempo es una sucesión de variables aleatorias ordenadas de acuerdo a una unidad de tiempo, \(Y_{1}, \ldots, Y_{T}\).
¿Por qué usar series de tiempo?
Definción:
\[ \Delta Y_{t-i} = Y_t-Y_{t-i} \]
Ejemplos:
\[ \Delta Y_{t} = Y_t-Y_{t-1} \]
Caso general:
\[ L^{j}Y_t=Y_{t-j} \]
Ejemplos:
\[ L^1Y_t=LY_t=Y_{t-1} \] \[ L^2Y_t=Y_{t-2} \]
\[ L^{-2}Y_t=Y_{t+2} \]
\[
L^{i}L^{j}=L^{i+j}=Y_{t-(i+j)}
\] # Manipulando ts en `R
IPCEcuador.csvuu <- "https://raw.githubusercontent.com/vmoprojs/DataLectures/master/IPCEcuador.csv"
datos <- read.csv(url(uu),header=T,dec=".",sep=",")
IPC <- ts(datos$IPC,start=c(2005,1),freq=12)
plot(IPC)
La serie tiene tendencia creciente. Tratemos de quitar esa tendencia:
plot(diff(IPC)) # Se puede ver una inlfacion estable
abline(h=0)
Se ha estabilizado, pero podemos hacerlo aún más con el logartimo de la diferencia:
plot(diff(log(IPC))) #Tasa de variacion del IPC
La serie no tiene tendecia y es estable. Ahora, si deseo trabajar con un subconjunto de datos, puedo…
# Solo quiero trabajar con los datos de agosto 2008
IPC2 <- window(IPC,start=c(2008,8),freq=1)
plot(IPC2)
# IPC de todos los diciembres
IPC.dic <- window(IPC,start=c(2005,12),freq=T)
plot(IPC.dic)
points(IPC.dic)
Si tengo mensuales y necsito trabajar con el IPC anual:
aggregate(IPC)
## Time Series:
## Start = 2005
## End = 2012
## Frequency = 1
## [1] 1213.09 1241.75 1262.13 1349.24 1399.26 1435.96 1491.87 1542.53
A continuación algunas transformaciones frecuentes y su interpretación:
| Transformación | Interpretación |
|---|---|
| \(z_t=\nabla y_t=y_t-y_{t-1}\) | Cambio en \(y_t\). Es un indicador de crecimiento absoluto. |
| \(z_t=ln(y_t)-ln(y_{t-1})\approx \frac{y_t-y_{t-1}}{y_{t-1}}\) | Es la tasa logarítmica de variación de una variable. Es un indicador de crecimiento relativo. Si se multiplica por 100 es la tasa de crecimiento porcentual de la variable |
| \(z_t=\nabla[ln(y_t)-ln(y_{t-1})]\) | Es el cambio en la tasa logarítmica de variación de una variable. Es un indicador de la aceleración de la tasa de crecimiento relativo de una variable. |
Veamos un gráfico más interesante usando un conjunto de datos anterior, vamos a:
estadísticas Turismo.csvts y graficaruu <- "https://raw.githubusercontent.com/vmoprojs/DataLectures/master/estadisticas%20Turismo.csv"
datos<-read.csv(url(uu),header=T,dec=".",sep=";")
attach(datos)
# Visitas a Areas Naturales Protegidas
# Sumar por mes y año
mensual<-aggregate(TOTALMENSUAL,by=list(mesnum,Year),FUN="sum") # Los datos sin mes es el total de ese anio
TOTALmensual<-ts(mensual[,3],start=c(2006,1),freq=12)
plot(TOTALmensual)
Se ve una tendencia creciente y también una cierta estacionalidad. Veamos la misma serie en gráficos más atractivos:
library(latticeExtra)
library(RColorBrewer)
library(lattice)
xyplot(TOTALmensual)
asTheEconomist(xyplot(TOTALmensual,
main="TOTAL VISITAS MENUSALES \n AREAS PROTEGIDAS")
,xlab="Year_mes",ylab="Visitantes")
Componentes
\[ Y_t = f(S_t,T_t,E_t) \] donde \(Y_t\) es la serie observada, \(S_t\) es el componente estacional, \(T_t\) es la tendencia y \(E_t\) es el término de error.
La forma de \(f\) en la ecuación anterior determina tipos de descomposiciones:
| Descomposición | Expresión |
|---|---|
| Aditiva | \(Y_t = S_t +T_t+E_t\) |
| Multiplicativa | \(Y_t = S_t *T_t*E_t\) |
| Transformación logaritmica | \(log(Y_t) = log(S_t) +log(T_t)+log(E_t)\) |
| Ajuste estacional | \(Y_t - S_t =T_t+E_t\) |
visitas.descompuesta<-decompose(TOTALmensual, type="additive")
plot(visitas.descompuesta)
Dentro de visitas.decompuesta tenemos los siguientes elementos:
$x = serie original$seasonal = componente estacional de los datos EJ: en marzo hay un decremento de 2502 (para cada dato)$trend = tendencia$random = visitas no explicadas por la tendencia o la estacionalidad$figure = estacionalidad (mismo que seasonal pero sin repetición)Visualmenete:
Forma alternativa de elegir: ver cuál es la que tiene un componente aleatorio menor.
Los datos en el archivo wine.dat son ventas mensuales de vino australiano por categoría, en miles de litros, desde enero de 1980 hasta julio de 1995. Las categorías son blanco fortificado (fortw), blanco seco (dryw), blanco dulce (sweetw), rojo (red), rosa (rose) y espumoso (spark).
direccion<- "https://raw.githubusercontent.com/dallascard/Introductory_Time_Series_with_R_datasets/master/wine.dat"
wine<-read.csv(direccion,header=T,sep="")
attach(wine)
head(wine)
## winet fortw dryw sweetw red rose spark
## 1 1 2585 1954 85 464 112 1686
## 2 2 3368 2302 89 675 118 1591
## 3 3 3210 3054 109 703 129 2304
## 4 4 3111 2414 95 887 99 1712
## 5 5 3756 2226 91 1139 116 1471
## 6 6 4216 2725 95 1077 168 1377
dulce <- ts(sweetw,start=c(1980,12), freq=12)
plot(dulce)
Tratemos la serie como un caso aditivo:
En funcion del grafico de la variable, se decide el “type” de la descomposicion La estacionalidad tiene valores negativos porque se plantea respecto de la tendencia
dulce.descompuesta<-decompose(dulce, type="additive")
plot(dulce.descompuesta)
a<-dulce.descompuesta$trend[27] # La tendencia era de 130 en mayo del 82
b<-dulce.descompuesta$seasonal[27] # El componente de la estacionalidad era este
c<-dulce.descompuesta$random[27] # Es el componente aleatorio
a+b+c # La sumatoria de la descomposicion de la serie da el valor real, si es aditiva
## [1] 127
Veamos el caso multiplicativo:
# Multiplicativa
dulce.descompuesta1<-decompose(dulce, type="multiplicative")
plot(dulce.descompuesta1)
a<-dulce.descompuesta1$trend[27] # La tendencia era de 130 en mayo del 82
b<-dulce.descompuesta1$seasonal[27] # El componente de la estacionalidad era este
c<-dulce.descompuesta1$random[27] # Es el componente aleatorio
a*b*c # La sumatoria de la descomposicion de la serie da el valor real, si es aditiva
## [1] 127
Veamos la forma alternativa de elección:
u1<-var(dulce.descompuesta$random,na.rm=T)
u2<-var(dulce.descompuesta1$random,na.rm=T)
cbind(u1,u2)
## u1 u2
## [1,] 1970.235 0.02247602
Se escoge la multiplicativa en este caso.
Las series se ofrecen generalmente sin tendencia ni estacionalidad. Veamos la serie sin tendencia:
dulce.detrended <- dulce-dulce.descompuesta$trend
plot(dulce.detrended)
abline(h=0)
Parece ser que hay un cambio en la varianza desde el 85.
Si descomponemos multiplicativamente en vez de restar se debe dividir.
plot(dulce-dulce.descompuesta$trend)
plot(dulce/dulce.descompuesta1$trend)
Existen formas de descomponer más sofisticadas, por ejemplo, usando la función stl.
dulce.stl<-stl(dulce,s.window="per")
plot(dulce.stl)
En este caso el calculo de la tendencia cambia, se calcula con formas no paramétricas. La barra del final es la desviacion estándar.
El método se resume en las fórmulas siguientes:
\[\begin{eqnarray} a_{t} &=& \alpha(x_{t}-s_{t-p})+(1-\alpha)(a_{t-1}+b_{t-1}) \nonumber \\ b_{t} &=& \beta(a_{t}-a_{t-1})+(1-\beta)b_{t-1} \nonumber \\ s_{t} &=& \gamma(x_{t}-a_{t})+(1-\gamma)s_{t-p} \nonumber \end{eqnarray}\]El método de Holt-Winters generaliza el método de suavizamiento exponencial.
Veamos un modelo más sencillo:
\[\begin{eqnarray} x_{t} &=& \mu_{t}+w_{t} \nonumber \\ \mu_{t}=a_{t}&=& \alpha x_{t} + (1-\alpha) a_{t-1} \nonumber \end{eqnarray}\]dulce.se <- HoltWinters(dulce,beta=0,gamma=0)
plot(dulce.se)
Es un suavizamiento HW sin tendencia y sin componente estacional. La serie roja son los datos con suavizamiento exponencial y la negra son los observados. R buscó el alpha que le pareció apropiado.
Usemos un alpha deliberado:
dulce.se1 <- HoltWinters(dulce,alpha=0.8,beta=0,gamma=0)
plot(dulce.se1)
¿Qué pasó con los errores?
dulce.se$SSE # Suma de los residuos al cuadrado (de un paso)
## [1] 963408.2
dulce.se1$SSE # Suma de los residuos al cuadrado (de un paso)
## [1] 1132577
Es decir, el criterio para la busqueda de los parámetros es la minimización del SSE.
Total mensual de pasajeros (en miles) de líneas aéreas internacionales, de 1949 a 1960.
data(AirPassengers)
str(AirPassengers)
## Time-Series [1:144] from 1949 to 1961: 112 118 132 129 121 135 148 148 136 119 ...
plot(AirPassengers)
Se aprecia tendencia y variabilidad. Podemos usar HW para predicción:
ap.hw<- HoltWinters(AirPassengers,seasonal="mult")
plot(ap.hw)
ap.prediccion <- predict(ap.hw,n.ahead=48)
ts.plot(AirPassengers,ap.prediccion,lty=1:2,
col=c("blue","red"))
n <-200
mu <- 0
sdt <- 3
w <- rnorm(n,mu,sdt)
¿Cómo se si algo tiene ruido blanco? : Analizo la función de autocorrelación muestral.
\[ \hat\rho_k = \frac{\hat\gamma_k}{\hat\gamma_0} \] donde \[ \hat\gamma_k = \frac{\sum(Y_t-\bar{Y})(Y_{t+k}-\bar{Y})}{n} \] \[ \hat\gamma_0 = \frac{\sum(Y_t-\bar{Y})^2}{n} \]
Se asume que \(\hat\rho_k \sim N(0,1/n)\)
acf(w)
Si se sale de las franjas, si hay correlación y no hay ruido blanco
Proceso estocástico autoregresivo de primer orden
\[\begin{eqnarray} (Y_{t}-\delta) = \alpha_{1}(Y_{t-1}-\delta)+u_{t} \nonumber \end{eqnarray}\]Donde \(\delta\) es la media de \(Y\) y \(u_{i}\) es un ruido blanco no correlacionado.
Trabajaremos con datos de \(M1\) (WCURRNS dinero en circuación fuera de los Estados Unidos) semanales de los Estados Unidos desde enero de 1975.
uu <- "https://raw.githubusercontent.com/vmoprojs/DataLectures/master/WCURRNS.csv"
datos <- read.csv(url(uu),header=T,sep=";")
names(datos)
## [1] "DATE" "VALUE"
attach(datos)
value.ts <- ts(VALUE, start=c(1975,1),freq=54)
ts.plot(value.ts)
Estacionariedad: La serie es estacionaria si la varianza no cambia
acf(value.ts)
Esta es la marca de una serie que NO es estacionaria, dado que la autocorrelación decrece muy lentamente.
plot(stl(value.ts,s.window="per"))
Una forma de trabajar con una serie esacionaria es quitarle el trend
valuetrend<- value.ts- stl(value.ts,s.window="per")$time.series[,2]
plot(valuetrend)
Reminder es lo que queda sin tendencia ni estacionalidad
valuereminder<-
stl(value.ts,s.window="per")$time.series[,3]
Veamos cómo quedo la serie:
acf(valuereminder)
Se puede decir que la hicimos una serie estacionaria
Otra forma de hacer estacionaria una serie es trabajar con las diferencias
acf(diff(value.ts))
Nos indica que hay una estructura en la serie que no es ruido blanco pero SI estacionaria (cae de 1 a “casi”" cero)
El siguiente paso es modelar esta estructura. Un modelo para ello es un modelo autorregresivo. Simular un \(AR(1)\).
y <- arima.sim(list(ar=c(0.99),sd=1),n=200)
plot(y)
acf(y)
¿Cuáles son los parámetros del arima.sim? Hemos simulado \(Y_t = \phi_{0}+\phi_{1}Y_{t-1} = \phi_{0}+0.99Y_{t-1}\).
Simulemos el modelo: \(Y_t = 0.5Y_{t-1} - 0.7Y_{t-2} + 0.6Y_{t-3}\)
ar3 <- arima.sim(n=200,list(ar=c(0.5,-0.7,0.6)),sd=5)
ar3.ts = ts(ar3)
plot(ar3.ts)
acf(ar3)
Las autocorrelaciones decaen exponencialemente a cero
Autocorrelaciones parciales: nos ayunda a determinar el orden del modelo.
La autocorrelación parcial es la correlación entre \(Y_t\) y \(Y_{t-k}\) después de eliminar el efecto de las Y intermedias.
En los datos de series de tiempo, una gran proporción de la correlación entre \(Y_t\) y \(Y_{t-k}\) puede deberse a sus correlaciones con los rezagos intermedios \(Y_1,Y_2,\ldots,Y_{t-k+1}\). La correlación parcial elimina la influencia de estas variables intermedias.
pacf(ar3)
ar(ar3)$aic
## 0 1 2 3 4 5
## 163.681281 162.849694 56.662316 0.000000 1.821970 2.933560
## 6 7 8 9 10 11
## 4.629183 6.515644 7.704677 9.682188 11.648622 12.895600
## 12 13 14 15 16 17
## 14.391168 13.820020 15.722147 17.717711 19.717700 21.354858
## 18 19 20 21 22 23
## 22.719834 24.497232 24.758384 26.601354 26.150993 28.107467
La tercera autocorrelación es la que esta fuera de las bandas, esto indica que el modelo es un AR(3)
Datos: precio de huevos desde 1901
uu <- "https://raw.githubusercontent.com/vmoprojs/DataLectures/master/PrecioHuevos.csv"
datos <- read.csv(url(uu),header=T,sep=";")
ts.precio <- ts(datos$precio,start=1901)
plot(ts.precio)
Veamos las autocorrelaciones:
par(mfrow=c(2,1))
acf(ts.precio)
pacf(ts.precio)
par(mfrow=c(1,1))
Las auto si decaen, no lo hacen tan rápido. No se puede decir si es estacionario o no.
Evaluemos un modelo:
modelo1 <- arima(ts.precio, order=c(1,0,0))
print(modelo1)
##
## Call:
## arima(x = ts.precio, order = c(1, 0, 0))
##
## Coefficients:
## ar1 intercept
## 0.9517 195.5066
## s.e. 0.0310 48.1190
##
## sigma^2 estimated as 712.3: log likelihood = -443.28, aic = 892.56
modelo1$var.coef
## ar1 intercept
## ar1 0.0009588394 -0.1532017
## intercept -0.1532017342 2315.4371739
¿Qué nos recomienda R?
ar.precio <- ar(ts.precio)
ar.precio
##
## Call:
## ar(x = ts.precio)
##
## Coefficients:
## 1
## 0.9237
##
## Order selected 1 sigma^2 estimated as 975.9
Analicemos los residuos
residuos = ar.precio$resid
# Los residuos debe estar sin ninguna estructura
par(mfrow = c(2,1))
plot(residuos)
abline(h=0,col="red")
abline(v=1970,col="blue")
acf(residuos,na.action=na.pass)
La prueba de Ljung-Box se puede definir de la siguiente manera.
\(H_0\): Los datos se distribuyen de forma independiente (es decir, las correlaciones en la población de la que se toma la muestra son 0, de modo que cualquier correlación observada en los datos es el resultado de la aleatoriedad del proceso de muestreo).
\(H_a\): Los datos no se distribuyen de forma independiente.
La estadística de prueba es:
\[ Q=n\left(n+2\right)\sum _{k=1}^{h}{\frac {{\hat {\rho }}_{k}^{2}}{n-k}} \]
donde n es el tamaño de la muestra, \(\hat\rho_{k}\) es la autocorrelación de la muestra en el retraso k y h es el número de retardos que se están probando. Por nivel de significación \(\alpha\), la región crítica para el rechazo de la hipótesis de aleatoriedad es \[ Q>\chi _{1-\alpha ,h}^{2} \]
donde \(\chi _{1-\alpha ,h}^{2}\) es la \(\alpha\)-cuantil de la distribución chi-cuadrado con \(m\) grados de libertad.
La prueba de Ljung-Box se utiliza comúnmente en autorregresivo integrado de media móvil de modelado (ARIMA). Tenga en cuenta que se aplica a los residuos de un modelo ARIMA equipada, no en la serie original, y en tales aplicaciones, la hipótesis de hecho objeto del ensayo es que los residuos del modelo ARIMA no tienen autocorrelación. Al probar los residuales de un modelo ARIMA estimado, los grados de libertad deben ser ajustados para reflejar la estimación de parámetros. Por ejemplo, para un modelo ARIMAa (p,0,q), los grados de libertad se debe establecer en \(h-p-q\).
Box.test(residuos,lag=20,type="Ljung")
##
## Box-Ljung test
##
## data: residuos
## X-squared = 20.526, df = 20, p-value = 0.4255
¿Es ruido blanco?